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The detection of gravitational waves from astrophysical sources is probably one of the most 
keenly awaited events in the history of astrophysics. The paucity of gravitational wave sources and 
the relative difficulty in detecting such waves, as compared to those in the electromagnetic domain, 
necessitate the development of optimal data analysis techniques to detect the signal, as well as to 
extract the maximum possible information from the detected signals. Coalescing binary systems are 
one of the most promising sources of gravitational waves. This is due to the fact that such sources 
are easier to model and thus one can design detection strategies particularly tuned to such signals. A 
lot of attention has been devoted in the literature studying such techniques and most of the work has 
revolved around the Weiner filtering and the maximum likelihood estimators of the parameters of the 
binary system. We investigate such techniques with the aid of differential geometry which provides 
geometric insight into the problem. Such a formalism allows to explore the merits and demerits of a 
detection scheme independent of the parameters chosen to represent the waveform. The formalism 
also generalises the problem of choosing an optimal set of templates to detect a known waveform 
buried in noisy data. We stress the need for finding a set of convenient parameters for the waveform 
and show that even after the inclusion of the second-order post-Newtonian corrections, the waveform 
can essentially be detected by employing a one-dimensional lattice of templates. This would be very 
useful both for the purpose of carrying out the simulations as well as for the actual detection process. 
After setting up such a formalism we carry out a Monte Carlo simulation of the detection process 
for the initial LIGO/ VIRGO configuration for the first post-Newtonian corrected coalescing binary 
waveform. We compare the results of our simulations with the currently available estimates of the 
accuracies in the determination of the parameters and the probability distribution of the maximum 
likelihood estimators. Our results suggest that the covariance matrix underestimates, by over a factor 
of two, the actual errors in the estimation of parameters even when the signal-to-noise ratio is as high 
as 20. As only a tiny fraction of the events is expected to be detected with a signal-to-noise higher 
than this value, the covariance matrix is grossly inadequate to describe the errors in the measurement 
of the parameters of the waveform. It is found from our Monte Carlo simulations that the deviations 
from the covariance matrix are more in the case of the first post-Newtonian waveform than in the 
case of the Newtonian one. Inclusion of higher-order post-Newtonian corrections introduces new 
parameters that are correlated with those at the lower post-Newtonian waveform. Such correlations 
are expected to further increase the discrepancy of the covariance matrix results with those inferred 
from Monte Carlo simulations. Consequently, numerical simulations that take into account post- 
Newtonian corrections beyond the first post-Newtonian order are needed in order to get a clearer 
picture about the accuracy in the determination of parameters. We find that with the aid of the 
instant of coalescence the direction to the source can determined more accurately than with the 
time of arrival. 



I. INTRODUCTION 

Laser interferometric detectors of gravitational waves such as the LIGO § and VIRGO 
are expected to be operational by the turn of the century. Gravitational waves from 
coalescing binary systems of black holes and neutron stars are relatively 'clean' waveforms 
in the sense that they are easier to model and for this reason they are amongst the most 
important candidate sources for interferometric detectors. Binary systems are also valu- 
able sources of astrophysical information as one can probe the universe up to cosmological 
distances. For instance, statistical analysis of several binary coalescences enables the esti- 
mation of the Hubble constant to an accuracy better than 10% Events that produce 
high signal-to-noise ratio can be potentially used to observe such non-linear effects, such as 
gravitational wave tails, and to put general relativity into test in the strongly non-linear 



1 



regime g . Due to the weak coupling of gravitational radiation with matter the signal wave- 
form is in general very weak and will not stand above the detector noise. In addition to 
the on-going efforts to reduce the noise, and hence increase the sensitivity of the detector, 
a considerable amount of research activity has gone into the development of efficient and 
robust data analysis techniques to extract signals buried in very noisy data. For a recent 
review on gravitational waves from compact objects and their detection see Thorne [^0. 

Various data analysis schemes have been suggested for the detection of the 'chirp' wave- 
form from such systems [^-[l0|] . Among them the technique of Weiner filtering is the most 
promising [p^|-p^. Briefly, this technique involves correlating the detector output with a 
set of templates, each of which is tuned to detect the signal with a particular set of param- 
eters. This requires the signal to be known to a high level of accuracy. The fully general 
relativistic waveform from a coalescing binary system of stars is as yet unavailable. In the 
absence of such an exact solution, there have been efforts to find solutions perturbatively. 
Most of the work done in this area aims at computing the waveform correct to a high de- 
gree of accuracy so that theoretical templates computed based on them will obtain a large 
signal-to-noise ratio (SNR) when correlated with the detector output if the corresponding 
signal is present. In general, the number of parameters increases as we incorporate the 
higher order corrections. It is clear that the number of templates depends upon the number 
of signal parameters. As a consequence, the computing power for an on-line analysis will be 
greater for a larger number of parameters. In view of this restriction in computing power it 
is necessary to choose the templates in an optimal manner. This paper in part deals with 
this question. Investigations till now have been restricted to either choosing a finite subset 
of the signal space as templates or choosing templates from the 'Newtonian' or the 

'first post-Newtonian' family of waveforms [p"5|-p8[. We generalize this problem by using 
the language of differential geometry. We show that it is unnecessary to restrict oneself to 
templates that are matched exactly to any particular signal. 

Differential geometry has been used in statistics before (see and references therein) 
and the standard approach is to treat a set of parameterised probability distributions corre- 
sponding to a particular statistical model as a manifold. The parameters of the distribution 
serve as coordinates on this manifold. In statistical theory one frequently comes across the 
Fisher information matrix whose inverse gives a lower bound for the errors in the estima- 
tion of the parameters of a distribution. The Fisher information matrix turns out to be a 
very natural metric on the manifold of probability distributions and this metric can be used 
profitably in understanding the properties of a particular statistical model. Here in our pa- 
per we treat the set of coalescing binary signals corresponding to various parameters of the 
binary as a submanifold in the linear space of all detector outputs. We show in this paper 
that both the above mentioned manifolds are equivalent as far as their metrical properties 
are concerned. The geometric approach turns out to be useful not only in clarifying various 
aspects of signal analysis but also helps us to pose the question of optimal detection in a 
more general setting. 

Once a signal has been detected we can estimate the parameters of the binary. We 
assume that the parameters of the signal are the same as those of the template with which 
the maximum correlation is obtained. The errors involved in such an estimation has been 
worked out by several authors |^,^ 20 -2^, for the case of 'high' SNR and for the Newtonian 
and post-Newtonian waveforms using a single and a network of detectors. For the case of 
low SNRs one has to resort to numerical simulations. We have started a project to carry 
out exhaustive numerical simulations specifically designed to compute the errors in the 
estimation of parameters and covariances among them at various post-Newtonian orders, 
for circular and eccentric orbits, with and without spin effects and for different optical 
configurations of the interferometer. In this paper we report the results for the case of the 
initial LIGO configuration taking into account only the first post-Newtonian corrections and 
assuming circular orbits. Going beyond this needs tremendous amount of computing power 
which is just becoming available. 

The rest of the paper is organized as follows. In section || we describe the waveform from 
a coalescing binary system at various post-Newtonian orders. We introduce, following pS] , 
a set of parameters called 'chirp times'. These parameters are found to be very convenient 
when we carry out Monte Carlo simulations. It turns out that the covariance matrix is 
independent of these parameters and hence it is sufficient to carry out the simulations only 
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for a particular set of parameters. In section [II we develop a geometric interpretation of 
signal analysis. We begin by introducing a metric on the manifold from a scalar product, 
which comes naturally from the theory of matched filtering and then show that this metric 
is the same as the one used by Amari jl9|. Using the geometric approach we address 
the question of optimal filter placement and show that for the purpose of detection it is 
optimal to choose the templates outside the signal manifold. The covariance matrix of 
errors and covariances is shown to be the inverse of the metric on the manifold. In section 
IV we discuss the results of our simulations and compare the numerically obtained values 
and those suggested by the covariance matrix. We find substantial discrepancies in the 
predictions of the two methods. It is believed that the coalescing binary waveform shuts 
off abruptly at the onset of the plunge orbit. This has a major effect on the computations 
of the covariance matrix as well as on the Monte Carlo simulations. We discuss the effects 
of higher post-Newtonian corrections to the waveform. We also emphasisize the use of the 
instant of coalescence as a parameter in order to determine the direction to the source rather 
than the time of arrival [pol . Finally in section ^ we summarise our results and indicate 
future directions. 



II. COALESCING BINARY WAVEFORMS 



For the purpose of constructing templates for on-line detection, it is sufficient to work 
with the so called restricted post-Newtonian gravitational waveform. In this approximation 
the post-Newtonian corrections are incorporated only in the phase of the waveform, while 
ignoring corresponding corrections to the amplitude |30|] . Consequently, the restricted post- 
Newtonian waveforms only contain the dominant frequency equal to twice the orbital fre- 
quency of the binary computed up to the relevant order. In the restricted post-Newtonian 
approximation the gravitational waves from a binary system of stars, modeled as point 
masses orbiting about each other in a circular orbit, induce a strain h(t) at the detector 
given by 



/l(i) = A(^/(t))2/3cos[(p(t)], 



(2.1) 



where f{t) is the instantaneous gravitational wave frequency, the constant A involves the 
distance to the binary, its reduced and total mass, and the antenna pattern of the detector 
pO| , and the phase of the waveform ip{t) contains several pieces corresponding to different 
post-Newtonian contributions which can be schematically written as 



(2.2) 



Here ifo(t) is the dominant Newtonian part of the phase and (/?„ represents the nth order 
post-Newtonian correction to it. In the quadrupole approximation we include only the 
Newtonian part of the phase given by pOl 



ip{t) = (fait) 



167r/aTo 



fa 



-5/3' 



(2.3) 



where f(t) is the instantaneous Newtonian gravitational wave frequency given implicitly by 

-8/3' 



t-ta=To 



1 



(2.4) 



To is a constant having dimensions of time given by 



(2.5) 



and fa and $ are the instantaneous gravitational wave frequency and the phase of the signal, 
respectively, at t — ta- The time elapsed starting from an epoch when the gravitational wave 
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frequency is fa till the epoch when it becomes infinite will be referred to as the chirp time 
of the signal. In the quadrupole approximation tq is the chirp time. The Newtonian part 
of the phase is characterised by three parameters: (i) the time of arrival ta when the signal 
first becomes visible in the detector, (ii) the phase of the signal at the time of arrival and 
(iii) the chirp mass M — {fjpM'^y/^, where /i and Af are the reduced and the total mass 
of the binary, respectively. At this level of approximation two coalescing binary signals of 
the same chirp mass but of different sets of individual masses would be degenerate and thus 
exhibit exactly the same time evolution. This degeneracy is removed when post-Newtonian 
corrections are included. 

When post-Newtonian corrections are included the parameter space of waveforms acquires 
an extra dimension. In this paper we show that even when post-Newtonian corrections up to 
relative order c~^, where c is the velocity of light, are included in the phase of the waveform 
it is possible to make a judicious choice of the parameters so that the parameter space 
essentially remains only three dimensional as far as the detection problem is concerned. It 
should, however, be noted that the evolution of the waveform must be known to a reasonably 
high degree of accuracy and that further off-line analysis would be necessary to extract useful 
astrophysical information. 

With the inclusion of corrections up to second post-Newtonian order the phase of the 
waveform becomes Pll 



(fit) ^ ifioit) + ipi{t) + ipi.r^{t) + ip2{t) 



(2.6) 



where ^po{t) is given by (2.3) and the various post- Newtonian contributions are given by 
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(2.7) 
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and 



ip2{t) = 87r/aT2 
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Now f{t) is the instantaneous gravitational wave frequency correct up to second post- 
Newtonian order given implicitly by 
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In the above equations r's are constants having dimensions of time which depend only on 
the two masses of the stars and the lower frequency cutoff of the detector fa- The total 
chirp time now consists of four pieces: the Newtonian contribution tq is given by (2.5) and 
the various post-Newtonian contributions are 
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where 77 = fj,/M. The phase ( p.6[ ) contains the reduced mass fi in addition to the chirp 
mass M . Taking {A4 , 77) to be the post-Newtonian mass parameters the total mass and the 
reduced mass are given by M = A4r]~^^^, fi = Mif/^. Note that in the total chirp time r of 
the signal the 1.5 post-Newtonian contribution appears with a negative sign thus shortening 
the epoch of coalescence. 

With the inclusion of higher order post-Newtonian corrections a chirp template is char- 
acterised by a set of four parameters which we shall collectively denote by A^, /x = 1, . . . , 4. 
At the first post-Newtonian approximation instead of working with the parameters — 
{ta, M, 77} we can equivalently employ the set {ta, 'I', tq, n} for the purpose of 
constructing templates. This, as we shall see later, has some advantageous. However, at 
post-Newtonian orders beyond the first we do not have a unique set of chirp times to work 
with. 

The parameters ta and $ are kinematical that fix the origin of the measurement of time and 
phase, respectively, while the Newtonian and the post-Newtonian chirp times are dynamical 
parameters in the sense that they dictate the evolution of the phase and the amplitude of 
the signal. It may be mentioned at this stage that in most of the literature on this subject 
authors use the set of parameters {tc, <&c, -M, 77} where tc is the instant of coalescence 
and $(7 is the phase of the signal at the instant of coalescence. In terms of the chirp times 
we have introduced, tc is the sum of the total chirp time and the time of arrival and is 
a combination of the various chirp times and $ : 

tc = ta + To + n - Ti,5 + T2; (2.14) 
$C = $ + ^^ro + AirfaTi - hnfaTi.^ + %TTfaT2. (2.15) 
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In the stationary phase approximation the Fourier transform of the restricted first-post- 

for positive 
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(2.16) 



is a normalisation constant, A'', /i = 1, . . . , 6, represent the various post-Newtonian param- 
eters 



A^ - {ia,$,ro,ri,ri.5,r2}, (2.17) 
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(2.19) 

(2.20) 
(2.21) 
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For / < the Fourier transform is computed using the identity h{—f) = h*{f) obeyed by 
real functions h{t). In addition to the above mentioned parameters we shaU introduce an 
ampUtude parameter A in section HI. 



III. A GEOMETRIC APPROACH TO SIGNAL ANALYSIS 



In this section we apply the techniques of di fferen tial geometry to the problem of detecting 
weak signals embedded in noise. 



In section 



III A 



we introduce the concept of the signal 
manifold and elaborate on the relationship of our approach with that of Amari |ig| ]. Our 
discussion of th e vec tor space of all detector outputs is modeled after the discussion given in 
p2| . In section IIIE| we deal with the problem of choosing a set of filters for o n-line analysis 
which would optimise the task of detection of the signal. In section [II C we deal with 
the dimensionality of the chirp manifold when we incorporate higher order post-Newtonian 
corrections. It is found, that due to covariances between the parameters, it is possible to 
introduce an effective dimension which is less than the dimension of the manifold. This has 
very important implications for the detection problem. 



A. Signal manifold 



The output of a gravitational wave detector such as the LIGO, will comprise of data 
segments, each of duration T seconds, uniformly sampled with a sampling interval of A, 
giving the number of samples in a single data train to be = T/A. Each data train can be 
considered as a A^-tuple (x", x^, . . . , x^'-^) being the value of the output of the detector 
at time lA. The set of all such A^ tuples constitutes a A^-dimensional vector space V where 
the addition of two vectors is accomplished by the addition of corresponding time samples. 
For later convenience we allow each sample to take complex values. A natural basis for 
this vector space is the time basis = (5*„ where m and i are the vector and component 
indices respectively. Another basis which we shall use extensively is the Fourier basis which 
is related to the time basis by a unitary transformation U: 

(3.1) 
(3.2) 

All vectors in V are shown in boldface, and the Fourier basis vectors and components of 
vectors in the Fourier basis are highlighted with a 'tilde'. 

In the continuum case each data train can be expanded in a Fourier series and will contain 
a finite number of terms in the expansion, as the output will be band limited. The expansion 
is carried out over the exponential functions exp (27rimt/T) which are precisely the Fourier 
basis vectors defined above. Though the index m takes both positive and negative values 
corresponding to positive and negative frequencies it is both, possible and convenient to 
allow m to take only positive values js^. Thus the vector space V can be considered as 
being spanned by the N Fourier basis vectors implying immediately that the number of 
independent vectors in the time basis to be also A^. This is the content of the Nyquist 
theorem which states that it is sufficient to sample the data at a frequency which is twice 
as large as the bandwidth of a real valued signal, where the bandwidth refers to the range 
of positive frequencies over which the signal spectrum is non zero. This factor of two does 
not appear in the vector space picture as we allow in general for complex values for the 
components in the time basis. 

A gravitational wave signal from a coalescing binary system can be characterised by a set 
of parameters A = (A°, A"'^, . . . , A^"^) belonging to some open set of the p-dimensional real 
space RP. The set of such signals s{t; A) constitutes a p-dimensional manifold S which is 
embedded in the vector space V. The parameters of the binary act as coordinates on the 



~ _ Tjmn _ 



-y 



e„ exp 



-^J^e^exp 



n=0 
W-1 



l-nimn 



N 
2Trimr 



6 



manifold. The basic problem of signal analysis is thus to determine whether the detector 
output vector x is the sum of a signal vector and a noise vector, x = s + n, or just the 
noise vector, x = n, and furthermore to identify which particular signal vector, among all 
possible. One would also like to estimate the errors in such a measurement. 

In the absence of the signal the output will contain only noise drawn from a stochastic 
process which can be described by a probability distribution on the vector space V. The 
covariance matrix of the noise C*-' is defined as, 



Cy=n»7i*i, (3.3) 

where an * denotes complex conjugation and an overbar denotes an average over an ensem- 
ble. If the noise is assumed to be stationary and ergodic then there exists a noise correlation 
function K{t) such that Cy = K{\i — j|A). In the Fourier basis it can be shown that the 
components of the noise vector are statistically independent and the covariance matrix 
in the Four ier basis will contain only diagonal terms whose values will be strictly positive: 
Cii — n^n". This implies that the covariance matrix has strictly positive eigenvalues. The 
diagonal elements of this matrix Cu constitute the discrete representation of the power 
spectrum of the noise Snif)- 

We now discuss how the concept of matched filtering can be used to induce a metric on 
the signal manifold. The technique of matched filtering involves correlating the detector 
output with a bank of filters each of which is tuned to detect the gravitational wave from a 
binary system with a particular set of parameters. The output of the filter, with an impulse 
response q, is given in the discrete case as 

N-l 

= 7^ E =^"9*" i-^^^rnn/N] . (3.4) 

V ^ n=0 

The SNR (p) at the output is defined to be the mean of C(m) divided by the square root of 
its variance: 



^ (3.5) 



(rn) ^(rn) ) 



1/2- 



By maximising p we can obtain the expression for the matched filter q(m) matched to a 
particular signal s{t; A^) as 

~„ g"(A^)exp[27r»mn/jV] 

9(m)l^ ) - 1 l-J-Oj 

where p has been maximised at the m"* data point at the output and where p ^ 1,2, ... ,p 
where p is the number of parameters of the signal. We now introduce a scalar product in 
V. For any two vectors x and y, 

(x,s)=5:C-xV=E^^4^- (3-7) 

In terms of this scalar product, the output of the matched filter q, matched to a signal 
s(A^), can be written as, 

C(„)(A'') = ^(x,s). (3.8) 



N 

As Cii is strictly positive the scalar product defined is positive definite. The scalar product 
defined above on the vector space V can be used to define a norm on V which in turn 
can be used to induce a metric on the manifold. The norm of a vector x is defined as 
|jx|| = (x, x)^^^. The ratio for the matched filter can be calculated to give p = (s, s)^^^. The 
norm of the noise vector will be a random variable (n, n)^^^ with a mean of ^/N as can be 
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seen by writing the expression for the norm of the noise vector and subsequently taking an 
ensemble average. 

The distance between two points infinitesimally separated on S can be expressed as a 
quadratic form in the differences in the values of the parameters at the two points: 

g^.dX'^dy = \\s{X>^ + dA^) - s(A^)|| = ll^'^A^II (3-9) 

= /^.|^V™'^ (3-10) 

The components of the metric in the coordinate basis are seen to be the scalar products of 
the coordinate basis vectors of the manifold. 

Since the number of correlations we can perform on-line is finite, we cannot have a filter 
corresponding to every signal. A single filter though matched to a particular signal will 
also detect signals in a small neighbourhood of that signal but with a loss in the SNR. The 
metric on the manifold quantifies the drop in the correlation in a neighbourhood of the 
signal chosen. Taking the output vector to be x and two signal vectors s(A) and s(A + dA) 
and using Schwarz's inequality we have, 

(x, s(A + dA)) - (x, s(A)) = (x, s(A + dA) - s(A)) < ||x|l |ls(A + dA) - s(A)|l (3.11) 

= ||x||5^,dA^dA^ (3.12) 

As is apparent the drop in the correlation can be related to the metric distance on the 
manifold between the two signal vectors. 

We now discuss Amari's ||l^] work in the context of using differential geometry in statistics 
and elaborate on the relationship with the approach we have taken. The set of parametrised 
probability distributions corresponding to a statistical model constitute a manifold. The 
parameterized probability distributions in the context of signal analysis of gravitational 
waves from coalescing binaries are the ones which specify the probability that the output 
vector will lie in a certain region of the vector space V given that a signal s(i; A) exists in 
the output which we denote as p(x|s(i; A)) Since it is not our intention to develop Amari's 
approach any further we will be brief and will make all the mathematical assumptions 
such as infinite differentiability of functions, interchangeability of the differentiation and 
expectation value operators, etc. 

The set of probability distributions p(x|s(A)) where A e i?^, constitutes a manifold V 
of dimension p. At every point on this manifold we can construct a tangent space on 
which we can define the coordinate basis vectors as 9^ = . Any vector A in this tangent 
space can be written as a linear combination of these coordinate basis vectors. We now 
define p random variables — — gpr log(p(x|s(A))). It can easily be shown that 5)7 = 0. 
We assume that these p random variables are linearly independent. By taking all possible 
linear combinations of these random variables we can construct another linear space T^. 
Each vector B in can be written as B = B^a^j,. The two vector spaces T° and 
are isomorphic to each other, which can be shown explicitly by making the correspondence 
<Jp, ^ dfj,. The vector space has a natural inner product defined on it which is the 
covariance matrix of the p random variables a^. This scalar product can be carried over to 
using the correspondence stated above. The metric on the manifold can be defined by 
taking the scalar product of the coordinate basis vectors 

9tiu = (9^, dv) = o^loy (3.13) 

In statistical theory the above matrix g^i, is called the Fisher information matrix. We will 
also denote the Fisher matrix, as is conventional, as F^jy. It is clearly seen that orthogonality 
between vectors in the tangent space of the manifold is related to statistical independence 
of random variables in T^. 

If we take the case of Gaussian noise the metric defined above is identical to the one 
obtained on the signal manifold by matched filtering. Gaussian noise can be described by 
the distribution. 
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where in the last step we have used the diagonal property of the matrix C*-' which implies 
that C-^ = l/Cu- 

As the noise is additive p(x|s(A)) can be written as p(x — s(A)). Assuming Gaussian noise 
we can write the expressions for the random variables cr^ as, 



1\ d / d 

'x - s(A), X - s(A)) = ( n, ^s(A) ) , 



2 / dXf' 
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where in the last step we have used x = s(A) + n. The covariance matrix for the random 
variables ct^ can be calculated to give 
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(3.16) 



which is the same metric as defined over the signal manifold. Thus, both the manifolds S 
and V are identical with respect to their metrical properties. We will henceforth restrict 
our attention to the signal manifold S. 

For the purpose of our analysis we will choose a minimal set of parameters character- 
izing the gravitational wave signal from a coalescing binary. We consider only the first 
post-Newtonian corrections. In section ^ we have already introduced the four parameters 
A'^ = {ta,^,TQ,Ti}. We now introduce an additional parameter for the amplitude and 
call it A" = A. The signal can now be written as s(/; A) = Ah{f;ta,^,TQ,Ti), where, 
A = {A, ta,^,Tf),Ti}. Numerically the value of the parameter A will be the same as that of 
the SNR obtained for the matched filter. We can decompose the signal manifold into a man- 
ifold containing normalised chirp waveforms and a one-dimensional manifold corresponding 
to the parameter A. The normalised chirp manifold can therefore be parameterized by 
{ta, 4", To, Ti}. This parameterization is useful as the coordinate basis vector ^ will be 
orthogonal to all the other basis vectors as will be seen below. 

In order to compute the metric, and equivalently the Fisher information matrix, we use 
the continuum version of the scalar product as given in j2^ , except that we use the two sided 
power spectral density. This has the advantage of showing clearly the range of integration 
in the frequency space though we get the same result using the discrete version of the scalar 
product. Using the definition of the scalar product we get 



df dS{f;\)drs*if;X) 
SM) dX^ d\- 



(3.17) 



Recall that in the stationary phase approximation the Fourier transform of the coalescing 

binary waveform is given by h{f) — A/'/~''/^exp « ^(/) — ^Mf): where 

given by equations ( ^.1^ - 2.23 ), fi — 1, ... ,4 and A'^ — {ta, <&, tq, ti}. Note, in 
particular, that in the phase of the waveform the parameters occur linearly thus enabling a 
very concise expression for the components of g^^. The various partial derivatives are given 

by 
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fVM(/)S(/;A), 



(3.18) 



where we have intr oduc ed t/jq ~ —i/A. On substituting the above expressions for the partial 
derivatives in eq. (3.17) we get, 







Snif) 





df 



(3.19) 
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The above definition of the amplitude parameter A, as in Culter and Flannagan |2|], disjoins 
the amplitude of the waveform from the rest of the par amete rs. Since V'o is pure imaginary 
and "ipfi's are real, it is straightforward to see from eq. ( [3.19 ) that 



500 = 1 and gop = 0; /i = 1, . . . , 4. (3.20) 

The rest of the components are seen to be independent of all the parameters except 
A i.e. g^i, cx A?. As A is unity for the normalised manifold the metric on the normalised 
manifold is flat. This implies not only that the manifold is intrinsically flat (in the stationary 
phase approximation) but also that the coordinate system used is Cartesian. If instead of 
the chirp time tq we use the parameter M then the metric coefflcients will involve that 
parameter and the coordinate system will no longer remain Cartesian. 



B. Choice of filters 



We now use the formalism developed to tackle the issue of optimal filter placement. Till 
now, it has been thought necessary to use a finite subset of the set of chirp signals as 
templates for detection. We show that this is unduly restrictive. We suggest a procedure 
by which the detection process can be made more 'efficient' by moving the filters out of the 
manifold. It must be emphasised that the algorithm presented below is both, simplistic and 
quite adhoc and is not necessarily the best. Moreover, we have implemented the algorithm 
only for the Newtonian case where the computational requirements are not very heavy. 
The signal manifold corresponding to higher post-Newtonian corrections will be a higher 
dimensional manifold and here the computational requirements will be substantial. The 
choice of optimal filters which span the manifold will then be crucial. 

Detection of the coalescing binary signal involves computing the scalar product of the 
output of the detector with the signal vectors. Subsequently one would have to maximise 
the correlations over the parameters and the number so obtained would serve as the statistic 
on the basis of which we can decide whether a signal is present in the given data train. 
Geometrically, this maximization corresponds to minimizing the angle between the output 
vector and the vectors corresponding to the normalised signal manifold. Using the cosine 
formula, 

(^W^^) l|s(A)|p + ||xf -||x-s(A)|p 

cosiO) = - — ,,,, ,,,, ~ — - — ,,,, ,,,, , (0.21) 

l|x||||s(A)|| 2||x||||s(A)|| 

it is clear that as ||s(A)|| is unity for the vectors belonging to the normalised signal manifold, 
maximising the scalar product is equivalent to minimising ||x— s(A)|| which is the distance 
between the tip of the output vector and the manifold. 

Given the constraints of computational power one would be able to evaluate only a finite 
number of these scalar products, say n^, in a certain amount of time depending on the data 
train length and the padding factor. It is therefore necessary to be able to choose the np 
filters in such a manner that the detection probability is maximal. We will need efficient 
on-line data analysis for two reasons: (i) To isolate those data trains which have a high 
probability of containing a signal and (ii) to determine the parameters of the binary early 
on during the inspiral and to use them for dynamical recycling techniques. 

Due to the finiteness of the filter spacing the signal parameters will in general not corre- 
spond to any of the n p filters chosen and this will lead to a drop in the maximum possible 
correlation. Till now attention has been focussed on identifying a set of optimal set of filters 
which are a discrete subset of the manifold. If detection is the sole purpose then the differ- 
ential geometric picture suggests that confining the filter vectors to the signal manifold is 
an unnecessary restriction and in fact non optimal. Thus it is worthwhile to explore making 
a choice of filters outside the manifold. The filter vectors will thus belong to V but will 
not, in general, correspond to any signal. It is, of course, true that we are sacrificing on the 
maximum possible correlation obtainable (when the signal's parameters coincide with those 
of the filter). Thus the problem essentially is to select np filter vectors which optimize the 
detection the efficiency of which depends upon the properties of the manifold. 
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In general a single filter vector would have to pick up signals over a region of the manifold. 
The extent of this region is determined by fixing a threshold on the correlation between the 
filter and any signal in the region. We will denote this threshold by 77, where takes a value 
which is close to, but less than unity. The typical value suggested for 77 is ^ 0.8, jl^. For a 
given filter q and a threshold 77 the region on the manifold corresponding to the filter will 
be denoted as Sq{ri), where 3^(1]) C S. Geometrically this region is the intersection of an 
open ball in V of radius 2^/^(1 — 77)"'^/^ (using the distance defined by the scalar product), 
with center q, and the manifold. The np filters taken together would have to 'span' the 
manifold which means that the union of the regions covered by each filter would be the 
manifold itself i.e., IJq'^q(^) — ^- ^^'^ filter q lies on the manifold, then the correlation 
function Cq(A) = (s(A),q) will reach its maximum value of unity in Sq{ri) when q = s(A) 
and will fall off in all directions. This means that the signals in the region which are further 
away from the filter are less likely to be picked up as compared to those in the immediate 
neighbourhood of the filter q. 

We assume that a finite subset of the normalised signal manifold has been chosen to act 
as filters by some suitable algorithm which taken together span the manifold. The 
number of filters will be determined by the available computing power. Consider one of 
these filters q, the region corresponding to it for a threshold of 77, iSq(77), and an arbitrary 
normalised vector qo which belongs to V but not necessarily S. By correlating the vector 
qo with vectors in Sq{ri) we obtain the correlation function Cq^(A) = (s(A),qo). We select 
that qo to serve as a more optimal filter which maximises the average of this correlation 
function: 

(s(A),qo)^^= ^ J (s(A),q„)V5dPA= / 1 s(AKA , qA , (3.22) 

Sc(7j) s<,('?) ' 

where g = det [g^u] ■ In the last step above the integration and the scalar product operations 
in the above equation have been interchanged. Moreover, for the normalised chirp manifold 
the metric does not depend upon the parameters in the coordinate system we have chosen 
and therefore, g^^, is a constant and the factor y/g cancels. We now use Schwarz's inequality 
to maximise the average correlation to obtain, 

qo=AA J s(A)dPA, (3.23) 

where Af is just a normalisation constant. 

We implemented the above algorithm for filter placement for the case of Newtonian signals 
with certain modifications. The normalised chirp waveform consists of three parameters 
($, ta, T„). If we keep ta and t„ fixed then the tip of the signal vector traces out a circle as 
we vary $. As any circle lies on a plane we can express a signal vector as a linear sum of 
two vectors where the two vectors differ only in the phase parameter and we take this phase 
difference to be tt/2. Thus we need only two mutually orthogonal filters to span the phase 
parameter. The time of arrival parameter ta is also a 'convenient parameter' as by the use 
of fast Fourier transforms the correlations for arbitrary time of arrivals can be performed 
at one go. It is therefore not profitable for us to maximise the average correlation over the 
phase parameter $ and the time of arrival ta. 

In view of the above restrictions we modified the the filter placement algorithm. We 
consider the correlation function for the case when the filter vector is on the manifold. We 
define the 'line of curvature' to be the curve on the manifold along which the correlation 
function falls the least. Figure (|^) illustrates the correlation function plotted as a function 
of To along the line of curvature. It is seen from the contour diagram of the numerically 
computed correlation function that the line of curvature lies nearly on the submanifold 
ta + Tn = constant, of the normalised chirp manifold. We take two curves passing through 
the point q in the region Sq{ri): 

1. ta + Tn — constant, <I> = and 

2. ta + Tn = constant, = tt/2. 
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We obtain one filter for each of the two curves by evaluating eq. 3.23 where the domains of 
integration correspond to the segments of the curves defined. 

Having determined the two filters we again plot the correlation function along the line of 
curvature as a function of tq in figure (^). The region of the manifold selected corresponds 
to a range of 5.8 sees to 6.0 sees in the parameter tq. It is interesting to note that we could 
have translated these values keeping the difference same without affecting our results. It 
can be seen that the correlation has a minimum at the center. In order to get a flatter 
correlation curve we select a linear combination of the original filter and the one obtained 
by averaging with suitable weights attached to each filter. This performs reasonably well 
as shown by the thick curve in the figure. The importance of having a flatter correlation 
function lies in the fact all the signals in a region can be picked up with equal efficiency 
and the drop in the maximum possible correlation can be compensated for by lowering the 
threshold. The average correlation obtained for the optimal filter is only marginally better 
than that obtained for the filter placed on the manifold. 

In the discussion above, we had started with a fixed number of filters np on the manifold 
and obtained another set of n^? filters which performs marginally better than the former 
set. Equivalently we can try to increase the span of each filter retaining the same threshold 
but reducing the number of filters required. In Figure [l] we observe that the optimal filter 
chosen spans the entire region considered with a threshold greater than 0.9, whereas the 
filter on the manifold spans about half the region at the same threshold. This indicates that 
by moving the filters out of the manifold in the above manner it may be possible to reduce 
the number of filters by a factor of two or so. One must however, bear in mind that the 
bank of filters obtained in this way are not optimal. There is scope to improve the scheme 
further. Our analysis is indicative of this feature. 



C. Effective dimensionality of the parameter space of a second order post-Newtonian 

waveform 



It has already been shown that the first post-Newtonian waveform is essentially one- 
dimensional . We argue in this subsection that even the second post-Newtonian waveform 
is essentially one-dimensional and a one-dimensional lattice suffices to filter the waveform. 

A Newtonian waveform is characterised by a set of three parameters consisting of the 
time of arrival, the phase of the signal at the time of arrival and chirp mass (equivalently, 
Newtonian chirp time). In this case, for the purpose of detection, one essentially needs to 
employ a one-dimensional lattice of filters corresponding to the chirp mass, the time of arrival 
being taken care by the fast Fourier transform algorithm and the phase being determined 
using a two-dimensional basis of orthogonal templates. When post-Newtonian corrections 
are included in the phase of the waveform the number of parameters increases apparently 
implying that one needs to use a two-dimensional lattice of filters corresponding to, say, 
the chirp and reduced masses (equivalently the Newtonian and post-Newtonian chirp times) 
which in turn means that the number of templates through which the detector output needs 
to be filtered goes up by several orders of magnitude. One of us (BSS) has recently shown 
that for the purpose of detection it is sufficient to use a one-dimensional lattice of filters 
even after first post-Newtonian corrections are included in the phase of the waveform and 
the relevant parameter here is the sum of the Newtonian and post-Newtonian chirp times. 
What happens when corrections beyond the first post-Newtonian order are incorporated in 
the phase of the waveform? 

The coalescing binary waveform is now available up to second post-Newtonian order 
|l|,|§. Blanchet et al. argue that the phase correction due to the second order post- 
Newtonian (2PN) term induces an accumulated difference of 10 cycles in a total of 16000. 
Consequently, it is important to incorporate the 2PN terms in the templates. When the 
2PN terms are included it is useful to consider that the full waveform is parameterised by 
three additional parameters, corresponding to the chirp times at the 1, 1.5 and 2PN order 
(cf. Sec H), as compared to the Newtonian waveform. Of course, as far as the detection 
problem is concerned there is only one additional parameter since the chirp times are all 
functions of the two masses of the binary. However, for the purpose of testing general 
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relativity one can consider each of the chirp times to be independent of the rest |5|,^ . Our 
problem now is to find the dimensionality of the parameter space of a 2PN waveform. To 
this end we consider the ambiguity function C{\ , A) which is nothing but the correlation 
function of two normalised waveforms one of whose parameters (A) are varied by holding 
the parameters of the other fixed (A ) : 

C(A',A) = (g(A'),g(A)); (q(A'), g(A')) = (g(A ), (A)) = 1. (3.24) 

It is useful to think of A to be the parameters of a template and A to be that of a signal. 
With this interpretation the ambiguity function simply gives the span of a filter in the 
parameter space. 

The ambiguity function for the full waveform is a four-dimensional surface since there 
arc four independent parameters. To explore the effective dimensionality of the parameter 
we consider the set of parameters to be {ta, ^,1111,7712}, where mi and m2 are the two 
masses of the binary. We have shown the contours of the ambiguity function maximised 
over ta and <& (since these two parameters do not explicitly need a lattice of templates) 
in Fig. ^. The template at the centre of the plot corresponds to a binary waveform with 
mi = m2 = I.4M0 and the signal parameters are varied over the entire astrophysically 
interesting range of masses: mi, 1112 € [1.4, 10]M©. From these figures we find that the 
ambiguity function is almost a constant along a particular line in the mi-m2 plane. This 
means that a template at the centre of the grid spans a relatively large area of the parameter 
space by obtaining a correlation very close to unity for all signals whose masses lie on the 
curve along which the ambiguity function roughly remains a constant. It turns out that the 
equation of this curve is given by 

To + Ti — T1.5 + T2 ~ const. (3.25) 

Let us suppose we begin with a two-dimensional lattice of filters corresponding to a certain 
grid (albeit, nonuniform) laid in the mi-m2 plane. Several templates of this set will have 
their total chirp time the same. Now with the aid of just one template out of all those that 
have the same chirp time we can effectively span the region that is collectively spanned by 
all such filters. More precisely, we will not have an appreciable loss in the SNR in replacing 
all templates of a given total chirp time by one of them. Consequently, the signal manifold 
can be spanned by a one-dimensional lattice of templates. 

IV. ESTIMATION OF PARAMETERS 

In this section we discuss the accuracy at which the various parameters of a coalescing 
binary system of stars can be estimated. All our results are for a single interferometer 
of the initial LIGO-type which has a lower frequency cutoff at 40 Hz. At present it is 
beyond the computer resources available to us to carry out a simulation for the advanced 
LIGO. In the first part of this section we briefly review the well known results obtained 
for the variances and covariances in the estimation of parameters using analytical methods. 
Analytical methods assume that the SNR is sufficiently large (the so-called strong signal 
approximation) and implicitly use a continuum of the parameter space. In reality, however, 
these assumptions are not necessarily valid and hence it is essential to substantiate the 
results obtained using analytical means by performing numerical simulations. In the second 
part of this section we present an exhaustive discussion of the Monte Carlo simulations we 
have performed to compute the errors and covariances of different parameters. As we shall 
discuss below the computation of errors using the covariance matrix is erroneous even at an 
SNR of 10-20. Our estimation of 1-cr uncertainty in the various parameters, at low SNRs, 
is substantially larger than those computed using the covariance matrix. However, for high 
values of the SNR (> 25-30) Monte Carlo estimation agrees with the analytical results. 
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A. Covariance matrix 



In recent years a number of authors have addressed issues related to the variances expected 
in the parameter estimation ||5|, [l8| , p0| -p7| , In the standard method of computing the variances 
in the estimation of parameters one makes the assumption that the SNR is so large that 
with the aid of such an approximation one can first construct the Fisher information matrix 
r^j^ and then take its inverse to obtain the covariance matrix C^jy. In the strong signal 
approximation the Fisher information matrix and the covariance matrix are given by 

5.. = r,. = ^|^,^); c,. = r-V- (4.1) 

As we have seen before the Fisher information and consequently the covariance matrix 
is block diagonal and hence there is no cross-talk, implying vanishing of the covariances 
between the amplitude and the other parameters. Consequently, we need not construct, for 
the purpose of Weiner filtering, templates corresponding to different amplitudes. 

The judicious choice of parameters has also allowed a very elegant expression for the Fisher 
information matrix. It is particularly interesting to note that the information matrix does 
not depend on the values of the various parameters, except for the amplitude parameter, and 
hence is a constant as far as the parameter space of the normalised waveforms is considered. 

For the purpose of numerical simulations it is convenient to choose the set A'^ = 
{A,ta,^,To,Ti} where A is the amplitude parameter, ta and $ are the time of arrival 
of the signal and its phase at the time of arrival, and tq and ti are the Newtonian and 
the post-Newtonian coalescence times. For noise in realistic detectors, such as LIGO, the 
elements of the Fisher information matrix cannot be expressed in a closed form and, for the 
set of parameters employed, it is not useful to explicitly write down the covariance matrix in 
terms of the various integrals since the errors and covariances do not have any dependence 
on the parameters. We thus evaluate the information matrix numerically and then take its 
inverse to obtain the covariance matrix. Instead of dealing with the covariance matrix C it 
more instructive to work with the matrix of standard deviations and correlation coefHcients 
D which is related to the former by 

D,. = ( .'^r^'Z" (4.2) 

where ct^ — D^^ is the l-cr uncertainty in the parameter A^. The off-diagonal elements of D 
take on values in the range [—1,1] indicating how two different parameters are correlated: 
For [i ^ DfjLv — 1, indicates that the two are perfectly correlated, = — 1 means that 
they are perfectly anticorrelated and D^i, = implies that they are uncorrelated. Since 
the information matrix is block-diagonal, the amplitude parameter is totally uncorrelated 
with the rest and thus an error in the measurement of A will not reflect itself as an error in 
the estimation of the other parameters and vice versa. In contrast, as we shall see below, 
Newtonian chirp time is strongly anticorrelated to post-Newtonian chirp time, which implies 
that if in a given experiment tq happens to be estimated larger than its true value then it is 
more likely that t\ will be estimated to be lower than its actual value. Such correlations are 
useful as far as detection is concerned since they tend to reduce the number of templates 
needed in filtering a given signal. On the other hand, strong correlations increase the volume 
in the parameter space to which an event can be associated at a given confidence level. It 
seems to be in general true that a given set of parameters do not satisfy the twin properties 
of having small covariances and reducing the effective dimension of the manifold for the 
purpose of filtering. We elaborate on this point below. 

Given a region in a parameter space, it is useful to know the proper volume (as defined 
by the metric) of the manifold corresponding to the said region. In choosing a discrete set 
of filters for the detection problem one has to decide upon the maximum allowable drop in 
the correlation due to the finite spacing. Once this is fixed, the number of filters can be 
determined from the total volume of the manifold. For the detection problem it is beneficial 
to have a small volume whereas if the waveform is parameterized in such a way such that the 
manifold corresponding to it covers a large volume, then one can determine the parameters 
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to a greater accuracy. As a simple example let us consider a two-dimensional toy model: 
A = {Ai,A2}. We compare different signal manifolds each corresponding to a different 
parameterizations of the waveform. We assume that the covariance matrix and its inverse, 
the Fisher information matrix, to be: 



/(Til cri2 
\0'l2 cr22 



and Taiy = 
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The volume of the manifold corresponding to a region K, of the parameter space is given as, 

722[l-e](iAidA2, (4.4) 
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where e = 712/(711722) is the correlation coefficient. It can be clearly seen that if the corre- 
lation coefficient is small, keeping the variances in the parameters, fixed then the volume of 
the manifold is maximal. Since the parameters tq and ti are highly anticorrelated the proper 
volume corresponding to the region reduces to zero showing that the effective dimensionality 
of the manifold is less. 

Though, in principle, the variances and covariances are independent of the chirp time, in 
reality there arises an indirect dependence since one terminates a template at a frequency 
/ = l/(6^/^7rM) (where M is the total mass of the binary) corresponding to the plunge 
radius at a = 6M. Therefore, larger mass binaries are tracked over a smaller bandwidth so 
much so that there is less frequency band to distinguish between two chirps of large, but 
different, total mass. Consequently, at a given SNR the error in the estimation of chirp 
time s is l arger for greater mass binaries. This is reflected by the fact that the integrals in 
eq. ( ^.17 ) are somewhat sensitive to the value of the upper cutoff. (This also explains why 
the errors in the estimation of the chirp and reduced masses are larger for greater mass 
binaries ||2^,|2^.) In the following we assume that the noise power spectral density is that 
corresponding to the initial LIGO for which a fit has been provided by Finn and Chcrnoff 
pl|. For an SNR of 10 the matrix D is given by 



/ 1.0 
8.37 0.999 
3.16 
\ 




(4.5) 



for the Newtonian signal, and by 

/l.O 




20.4 0.997 
6.7 





-0.972 
-0.954 
45.1 



V 



0.911 ' 
0.881 
-0.982 
25.98 / 



(4.6) 



for the first post-Newtoni an co rrected signal. While computing varianaces and covariances 
the integrals in equation ( 3.19| ) are evaluated by chosing a finite upper limit of 1 kHz. In 
the above matrices the entries are arranged in the order {^, ia,$,ro} in the Newtonian 
case, {A,ta,^,To,Ti} in the post-Newtonian case, off diagonal elements are dimensionless 
correlation coefficients and, where appropriate, diagonal elements are in ms. The values 
quoted in the case of the Newtonian waveform are consistent with those obtained using a 
different set of parameters by Finn and Chernoff (2l|. In Fig. ^ we have plotted cr's, at an 
SNR of 10, as a function of the upper frequency cutoff fc for Newtonian and post-Newtonian 
chirp times and the instant of coalescence tc- We see that a is larger for higher mass binaries 
but this is because we have fixed the SNR. However, if we consider binaries of different total 
masses, all located at the same distance, then a more massive binary produces a stronger 
SNR so that in reality it may be possible to determine its parameters more accurately than 
that of a lighter binary. In Fig. |j we have plotted tr's for binaries all located at the same 
distance as a function of total mass. We fix one of the masses at a value of 1.4 Mq and 
vary the other from 1.4 Mq to 10 Mq. In computing the ct's plotted in this figure we have 
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terminated the waveform at the plunge orbit and normalised the SNR of a lOAf0-1.4M0 
binary system to 10. As a function of M the uncertainties in tq and ti initially falloff since 
the increase in the SNR for larger mass binaries more than compensates for the drop in the 
upper frequency cutoff. However, for M larger than a certain Mq the increase in SNR is 
not good enough to compensate for the drop in /c so much so that the uncertainties in tq 
and Ti increase beyond Mq. The parameter ta, however, falls off monotonically. 

B. Monte Carlo estimation of parameters 

In this Section we present the first in a series of efforts to compute the covariance matrix 
of errors through numerical simulations for a coalescing binary waveform at various post- 
Newtonian orders. Analytical computation of the covariance matrix, as in the previous 
section, gives us an idea of the covariances and variances but, as we shall see in this Section, 
at low SNR's it grossly underestimates the errors. Quite apart from the fact that the 
assumptions made in deriving the covariance matrix might be invalid at low SNR's, in a 
realistic detection and data analysis, other problems, such as discreteness of the lattice 
of templates, finite sampling of the data, etc., do occur. It therefore seems necessary to 
check the analytical calculations using numerical simulations to gain further insight into 
the accuracy at which physical parameters can be measured. This Section is divided into 
several parts: In the first part we highlight different aspects of the simulation, in the second 
part we briefly discuss the choice of templates for the simulation, in the third we elaborate 
on the Monte Carlo method that we have adopted to carry out our simulations and in the 
fourth we discuss problems that arise in a numerical simulation. The results of our study 
are discussed in the next Section. 



1. Parameters of the simulation 

Let s{t) be a signal of strength A characterised by a set of parameters A 

s{t-\)=Ah{t]\), {h,h)^l. (4.7) 

In data analysis problems one considers a discrete version {s*^|fc = 0, . . . ,N — 1} of the 
waveform s{t) sampled at uniform intervals in t : 

s^ = s(fcA); k = 0,...,N-l, (4.8) 

where A denotes the constant interval between consecutive samples and N is the total 
number of samples. The sampled output x'^ of the detector consists of the samples of the 
noise plus the signal: 

x'' = n'' + s''. (4.9) 

The sampling rate, fg = (also referred to as the sampling frequency) is the number 

of samples per unit time interval. In a data analysis problem the sampling frequency is 
determined by the signal bandwidth. If B is the signal bandwidth, i.e., if the Fourier 
transform of the signal is only nonzero over a certain interval B, then it is sufhcient to 
sample at a rate fs — 2B. In our case there is a lower limit in the frequency response of the 
detector since the detector noise gets very large below a seismic cutoff at about 10-40 Hz. 
As mentioned in the last Section there is also an upper limit in frequency up to which a chirp 
signal is tracked since one does not accurately know the waveform beyond the last stable 
orbit of the binary. This corresponds to gravitational wave frequency fc = l/(6'^/^7rAf ). For 
a neutron star-neutron star (ns-ns) binary fc ^ 1525 Hz while for a ns-black hole (of 10 Mq) 
(ns-bh) binary fc ^ 375 Hz. Due to constraints arising out of limited computational power, 
we terminate waveforms at 750 Hz even when fc is larger than 1000 Hz. Such a shutoff is 
not expected to cause any spurious results since, even in the case of least massive binaries 
of ns-ns, which we consider in this study, more than 99 % of the 'energy' is extracted by the 
time the signal reaches 750 Hz. We have carried out simulations with two types of upper 
cutoff: 
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1. one in which all templates, irrespective of their total mass, are shutoff beyond 750 Hz., 



2. a second in which the upper frequency cutoff is chosen to be 750 Hz or /c, whichever 
is lower. 

Consistent with these cutoffs the sampling rate is always taken to be 2 kHz. (We have 
carried out simulations with higher sampling rates and found no particular advantage in 
doing so nor did we find appreciable changes in our results.) 

In all our simulations, as in the previous section, we take the detector noise power spectral 
density S to be that corresponding to initial LIGO |2^ . For the purpose of simulations we 
need to generate noise corresponding to such a power spectrum. This is achieved by the 
following three steps: 

1. generate Gaussian white noise n''' with zero mean and unit variance. 



where an overbar denotes average over an ensemble, 
2. compute its Fourier transform 



V n" exp(27 

^ 1=0 



ikl/N) 



and 

3. multiply the Fourier components by the square root of the power spectral density. 



The resultant random process has the requisite power spectrum. In the above, the second 
step can be eliminated since the Fourier transform of a Gaussian random process is again 
a Gaussian, but with a different variance. In other words we generate the noise directly in 
the Fourier domain. The simulated detector output, in the presence of a signal s*"', in the 
Fourier domain is given by 

f = n'' + (4.10) 
where s*^ is the discrete Fourier transform of the signal. 



2. Choice of templates 

To filter a Newtonian signal we employ the set of parameters {ta,^,To} and to filter a 
post-Newtonian signal we employ the set {to, $, tq, n}. Templates need not explicitly be 
constructed for the time of arrival since computation of the scalar product in the Fourier do- 
main (and the availability of fast Fourier transform (FFT) algorithms) takes care of the time 
of arrival in essentially one computation (A^log2 N operations as opposed to N'^ operations, 
where N is the number of data points). Moreover, there exists a two-dimensional basis for 
the phase parameter which allows the computation of the best correlation with the aid of 
just two filters. Consequently, the parameter space is essentially one-dimensional in the case 
of Newtonian signals and t wo-dimensional in the case of post- Newtonian signals. (However, 



as shown in Section HI C it is to be noted that for the purpose of detection the effective 
dimensionality of the parameter space, even with the inclusion of second post-Newtonian 
corrections, is only one-dimensional.) We adopt the method described in Sathyaprakash 



and Dhurandhar |13| to determine the templates needed for chirp times. As described in 



3,E8| filters uniformly spaced in tq and ri covers the parameters space efficiently. 
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3. Monte Carlo method 



In order to compute variances and covariances numerically, we employ the Monte Carlo 
method. The basic idea here is to mimic detection and estimation on a computer by per- 
forming a very large number of simulations so as to minimize the uncertainties induced by 
noise fluctuations. In our simulations we generate a number of detector outputs {x'^} each 
corresponding to a definite signal s'^(A) of a certain strength, but corresponding to different 
realisations of the random process {n*'}. Computation of the covariance matrix involves fil- 
tering each of these detector outputs through an a priori chosen set (or lattice) of templates 
{q{t] tAfc)|fc = 1, . . . , Uf}, where Uf denotes the number of templates. The templates of the 
lattice each has a distinct set of values of the test parameters tAfc and together they span a 
sufficiently large volume in the parameter space. The simulated detector output is correlated 
with each member of the lattice to obtain the corresponding filtered output C(tAfc) : 

C(tAfc) = (a;,g(tAfe)). (4.11) 

For a given realisation of noise a particular template obtains the largest correlation and its 
parameters are the measwrerf values mA of the signal parameters. Thus, the measured values 
of the parameters are defined by 

maxC(tAfc) =C(^A). (4.12) 
fe 

The measured values, being specific to a particular realisation of noise, are random variables. 
Their average provides an estimation eA of the true parameter values and their variance is 
a measure of the error ax in the estimation: 



,A = ^A, al={^X-^\)\ £>^- = "'^ , (M5^i/), (4.13) 

where are the correlation coefficients between parameters A^ and A,^. In order to accu- 
rately determine ax a large number of simulations would be needed. If the measured values 
mA obey Gaussian statistics then after trials the variance is determined to a relative ac- 
curacy l/\/]Vt and estimated values can differ from their true values by ax/VNt. We have 
performed in excess of 5000 trials, for each input signal, and thus our results are accurate to 
better than 1 part in 70. Even more crucial than the number of simulations is the number of 
templates used and their range in the parameter space. We discuss these and other related 
issues next. 

The actual templates chosen, say for the parameter tq, in a given 'experiment' depend 
on the true parameters of the signal, the number of noise realisations employed and the 
expected value of the error. Let us suppose we have a first guess of the error in to, say arg- 
Then, we choose 51 uniformly spaced filters around fo (where tq is the signal chirp time) 
such that: 

tro e [fo - 5aro , tq + 5aro] ■ (4-14) 

This implies that we are covering a 5a width in tq at a resolution of cr^/S. The probability, 
that a template between 40- and 5a from the true signal 'clicks', being ^ 6 x 10^'', we 
are on safe grounds since, in a given simulation, we consider no more than 5000 trials. 
(In comparison, the probability that a template between 3a and 4<t clicks is 2.2 x 10~^ 
corresponding to an expected 13 events in 5000 trials.) For a post- Newtonian signal, which 
in effect needs to be spanned by a two-dimensional lattice of filters, the above choice of 
templates implies a requirement of 2601 x 2 filters in all. Here a factor of 2 arises because 
for each filter in the tq- ti space we will need two templates corresponding to the two 
independent values of the phase $: and 7r/2. In the case of a Newtonian signal, the lattice 
being one-dimensional, one can afford a much higher resolution and range. Even with the 
aid of a mere 201 templates we can probe at a a/10 resolution with a lOa range. 

We start off a simulation with the pretension that there is no knowledge of what the ax^s 
are. Thus, we choose as our first trial a very large ax and lay the lattice of templates. 
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With this lattice we perform a test run of 400 trials and examine the distribution of the 
measured values. If the distribution is not, as expected, a Gaussian then we alter ay. we 
decrease it if the distribution is too narrow and increase it if the distribution is too wide and 
does not show the expected falloff. In particular, we make sure that the templates at the 
boundary of the chosen range do not click even once and the skewness of the distribution is 
negligible. When for a certain ax a rough Gaussian distribution is observed then we carry 
out a simulation with a larger number of trials (typically 5000). We subject the measured 
values in this larger simulation to the very same tests described above. We only consider for 
further analysis such simulations which 'pass' the above tests and determine the estimates, 
variances and covariances of the parameters using the measured values, with the aid of 
equation (4.12). 



^. Numerical errors and remedies 



There are several sources of numerical errors that tend to bias the results of a simulation 
unless proper care is exercised to rectify them. In this Section we point out the most impor- 
tant ones and show how they can be taken care of. Due to memory restrictions, the present 
version of our codes work with single precision except the FFT, which is implemented in 
double precision. In future implementations we plan to carry out all computations in double 
precision. This will possibly reduce some of the numerical noise that occurs, especially at 
high SNRs, in the present simulations. 

1. Orthonormality of filters: For the sake of simplicity it is essential that the filters are 
normalised in the sense that their scalar product is equal to unity: {q,q) ~ 1- A 
waveform is normalized numerically using the discrete version of the scalar product: 

Af^ « 1 \ ? (4.15) 



As mentioned earlier we use a two-dimensional basis of filters for the phase parameter. 
Choosing the two filters to be orthogonal to each other makes the maximisation over 
phase easier. However, here care must be exercised. Two filters q{t;ta,TQ,Ti,^ = 0) 
and q{t] ta, tq, ti, $ = 7r/2) are apparently orthogonal to each other. The numerically 
computed 'angle' between the two filters, chosen in this manner, often turns out to be 
greater than ~ 10~^ radians. Consequently, one obtains erroneous correlations. In or- 
der to circumvent this problem we first generate two filters that are roughly orthogonal 
to each other, as above, and then use the Gram-Schmidt method to orthogonalise the 
two vectors. If an explicit numerical orthogonalisation such as this is not implemented 
then the measured values of the various parameters show spurious oscillations in their 
distribution and the estimated values of the parameters tend to get biased. 

2. Correlation function: The scalar product of two normalised templates q(i;tA) and 
g(i;tA') is given by 

C(t A^, t a:,) = {q{t; tA^), g(t; t A'J) , {q{t; tA), q{t; tA)) - (q(t; tA'), q{t; tA')) = 1, 

(4.16) 



where we have indicated the dependence of the scalar product on the various param- 
eters by explicitly writing down the parameter subscripts. Let us fix the parameters 
of one of the templates, say tA, and vary the parameters of the other template. Of 
particular interest is the behaviour of C maximised over all but one of the parameters, 
say A^o : 

G„ax(tA^,tA:,J = maxC(tA^,tA'J. (4.17) 

t 
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Ciaax is expected to drop monotonically as \Xuo — A^^| increases. However, we have 
observed departures from such a behaviour possibly arising out of numerical noise. 
Such a behaviour causes bias in the estimation of parameters, and consequently in the 
determination of their covariances, especially at high SNRs. We have found no remedy 
to this problem and some of our results at high SNRs may have biases introduced by 
this effect. (Sampling the templates at a higher rate did not help in curing this 
problem.) 

3. Grid effects: The parameters of a signal chosen for the purpose of simulation and 
detection can in principle be anything and in particular it need not correspond to any 
of the templates of the lattice. However, in practice we find that whenever the signal 
parameters do not correspond to a member of the lattice then the resultant simulation 
has a bimodal distribution of the measured values. This is, of course, expected since a 
signal not on the grid is picked up by two nearest templates along each direction in the 
parameter space. Sometimes we do find that the peaks corresponding to the bimodal 
distribution does not belong to the nearest neighbour filters but slightly away. This is 
related to the fact that the correlation function maximised over the time of arrival and 
the phase of the signal falls off much too slowly along the tq-ti direction and small 
deviation from a monotonic fall can cause biases. (Such biases would be present in the 
case when a signal corresponds to one of the grid points though the magnitude of the 
effect would be lower.) In order to avoid this problem, and the consequent shifts in the 
estimation of parameters and errors in the determination of variances and covariances, 
we always choose the parameters of the signal to be that corresponding to a template. 

4. Upper frequency cutoff and its effect on parameter estimation 

The Fisher information matrix computed using the stationary phase approximation in 



Section III does not include the effect of truncating the waveform at a = 6Af — the 
plunge cutoff. As mentioned before, we have carried out simulations for both with, 
and without, incorporating the upper cutoff. As the covariance matrix incorporating 
the upper cutoff is not available we have been able to compare the Monte Carlo results 
with the covariance matrix only for the latter case, where the cutoff is held fixed at 
750 Hz. If we incorporate the upper cutoff into the Monte Carlo simulations the errors 
in the parameters are reduced drastically. The effect of the upper cutoff is expected to 
be more important for the higher mass binaries such as the ones we have considered. 
The ambiguity function, in this case, no longer remains independent of the point on the 
manifold. In other words, the correlation between two chirps depends not only on the 
difference between the parameters of the signals, but also on the absolute values of the 
parameters. The correlation surface also ceases to be symmetric i.e. the correlation 
between two chirps also depends on the sign of (5A, where, 6X is the difference in 
the values of the parameters. As the computational power required for carrying out 
simulations for lower mass binaries is not available to us the simulations have been 
restricted to ns-bh star binaries, where the effect of the upper cutoff is important. 

5. Boundary effects 

For the purpose of simulations a grid of filters has to be set up 'around' the signal. 
The grid must be large enough so that the estimated parameters do not overshoot the 
boundary of the grid. This causes a problem as every value in the {tq, ti} plane does 
not lead to a meaningful value for the masses of the binary system. This does not 
however prevent us to construct a waveform with such a value for {tq, ti} even though 
the signal in general does not correspond to any 'real' binary system. This is valid, 
and even necessary, if we are to compare the numerical results with the covariance 
matrix. 

6. Incorporating the cutoff in the presence of boundary effects 

If we wish to incorporate the effects of the upper cutoff in simulations then we run into 
a serious problem, as we would have to know the total mass of the binary in order to 
compute the upper cutoff. For an arbitrary {to,ti} we can end up with negative and 
even complex values of the total mass and hence the upper cutoff at a = 6Af is not 
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meaningful. Thus, we cannot even construct a waveform for an arbitrary combination 
of {tqjTi}. Therefore, in such cases, we restrict ourselves to simulations where the 
grid lies entirely within valid limits for {tq, ti}. 



C. Results and Discussion 



Our primary o bjectiv e is to measure the variances and covariances following the method 
described in Se c. IV B 3 and study their departure from that predicted by analytical means 
(cf. Sec. [V A). We have carried out simulations for several values of the masses of the 



binary and in each case the signal strength (which is a measure of the SNR) is varied in the 
range 10-40. However, since the variances and covariances are independent of the absolute 
values of the parameters, for the parameter set th at we employ, results are only quoted 
corresponding to a typical binary system. (See Sec. IVA for a discussion of the covariance 
matrix.) Similar results are obtained in other cases too. We use two sets of parameters to 
describe our results. Monte Carlo simulations allow us to directly measure the amplitude, 
the time of arrival, the phase at the time of arrival and the chirp time(s). This is the set 

{A, ta,^,To,Ti}. 

As we shall see below, the instant of coalescence can be measured much more accurately 
than the time of arrival. As a consequence of this, the direction to the source can be 
determined at a much greater accuracy by employing tc as a parameter instead of ta- Thus, 
we also quote estimates and errors for the parameter tc- Since the error in the estimation 
of the phase is quite large, even at high SNRs, we ignore it in our discussions. 

We first deal with the Newtonian signal and highlight different aspects of the simulation 
and discuss the results at length. We then consider the first post-Newtonian corrected signal. 



1. Newtonian signal 

In the case of Newtonian signals the parameter space is effectively one-dimensional and, as 
mentioned earlier, in this case the lattice of templates covers a 10-cr range of the parameters 
at a resolution of cr/lO centred around the true parameters of the signal. 

In Fig. I we have shown the error ax in the estimation of parameters ta, tq and tc, as 
a function of SNR, deduced using the covariance matrix as solid lines and computed using 
Monte Carlo simulations as dotted lines. The curve corresponding to the covariance matrix 
is obtained using an upper frequency cutoff = 750 Hz consistent with that used in our 
simulations. The errorbars in the estimation of cta's are obtained using 4 simulations, each 
with 4000 trials. At low SNR's ax's have a larger uncertainty, as expected, and for p > 30 
this uncertainty is negligible, and sometimes smaller than the thickness of the curves, except 
in the case of at^ (see below, for a possible explanation). 

At low SNRs (10-15) there is a large departure of the various ua's from that inferred 
using the covariance matrix. At an SNR ~ 17 the two curves merge (except in the case 
of Ctc) indicating the validity of the covariance matrix results for this and higher SNR's. 
Interestingly, the agreement between Monte Carlo simulation results and those obtained 
using the covariance matrix is reached roughly at the same SNR irrespective of the parameter 
in question. 

We note that in spite of the fact that the time of arrival and the chirp time have large 
errors, the instant coalescence can be estimated very accurately — an order of magnitude 
better than either. What is puzzling, however, is that, in the case of tc, the Monte Carlo 
curve drops below the covariance matrix curve above an SNR of 15 and the two curves do 
not seem to converge to one another even at very high SNR's. Coincident with the crossover 
of the two curves, the error in the estimation of at^ increases, contrary to what happens for 
the other parameters, signalling that there is a large fluctuation in the estimation of at^- 
This behaviour, we guess, is an artifact of the low value of the sampling rate. Of course, 
our sampling rate is sufficiently high to respect the sampling theorem. However, since tc is 
determined to an accuracy an order of magnitude better than either ta or tq , a much higher 
resolution in template-spacing would be needed for determining the error in the instant 
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of coalescence than that used for estimating the errors in the time of arrival or the chirp 
time(s). Testing this claim, unfortunately, is beyond the computer resources at our disposal 
since we would need a sampling rate of about 10 kHz with a filter-spacing 10"^ s. We hope 
to be able to resolve this issue in course of time. Nevertheless, the fact that the error in the 
estimation of at^ first decreases with SNR and increases only after the two curves crossover, 
hints at the above possibility as a cause for this anomalous behaviour. This effect is also 
observed in the case of a post-Newtonian signal. 

In Table we have given the actual signal parameters A, estimated values of the parame- 
ters oA (cf. equation ( 4.13| )) and the corresponding errors in their estimation cta, for several 



values of the SNR. Errors inferred from the covariance matrix can be read off from Fig. |^. 
The estimated values are different from the true values, some of them being overestimated 
and some others being underestimated. However, the deviations are often larger than what 
we expect. In a simulation that uses Nt trials the estimated parameters e\ assuming a 
Gaussian distribution for the measured parameters mA, can be different from the true val- 
ues by (Tx/\/Nt. (In contrast, the measured values mA can differ from their true values by 
ax or more.) However, we often obtain a slightly larger deviation. 



2^< 



A A 



<34^, (4.18) 



and we are unable to resolve this discrepancy. A more concrete test for the simulations is 
the histogram n[^\) of the measured parameters, namely the frequency at which a given 
parameter clicks in a simulation. This is shown plotted in Fig. |^ for an SNR of 10. The 
skewness of the measured value is less than 10~^. These results lend further support to the 
Monte Carlo simulations. There are visible asymmetries in the distributions of tq and ta and 
the asymmetries in the two cases are of opposite sense. This can, of course, be understood 
from the fact that ta and tq have a negative correlation coefficient. The histogram of tc^ 
even at an SNR of 10, has very few non-zero bins. This of course reflects the fact that it 
is determined very accurately. We are unable to resolve the central peak in n{tc) since, as 
mentioned earlier, the sampling rate and resolution in tq are not good enough to do so. 



2. Post- Newtonian signal 

As opposed to the Newtonian case here we have essentially a two-dimensional lattice of 
filters corresponding to tq and ti . For the purpose estimating variances and covariances we 
lay a mesh consisting of 2601 x 2 un iformly spaced filters around the true parameters of the 



signal. As pointed out in Sec. HI B not all filters in the mesh, unlike in the Newtonian case, 
would correspond to the waveform from a realistic binary but that does not preclude their use 
in the Monte Carlo simulations. We shall see that the results of our simulations lend further 
support to the claim that for the purpose of detection, the parameter space need only be one- 
dimensional The results obtained for the first post-Newtonian signal are qualitatively 
similar to that of a Newtonian signal and we refer the reader, where appropriate, to the 
Newtonian case for a more complete discussion. 

In Fig. 1^ we have shown the error in the estimation of parameters tq, ti, tc and ta^ 
clockwise from top left, respectively, as a function of SNR. The solid and dotted curves are 
as in Fig. ||. Here again the upper frequency cutoff is taken to be 750 kHz. Just as in 
the case of a Newtonian signal here too the results obtained from Monte Carlo simulation 
are much higher than those obtained by employing the covariance matrix. At an SNR of 
10, which is the expected value for a majority of events that initial LIGO and VIRGO 
interferometers will observe, the Monte Carlo values are more than thrice as much as their 
corresponding covariance matrix values and at an SNR of 15 the errors are roughly twice 
that expected from the covariance matrix. In absolute terms, however, the errors are still 
quite small compared to the actual parameter values: for a ns-ns binary, at an SNR of 10, 

^-2.4%, ^^9.4%. (4.19) 

To Tl 

At an SNR of 10 the time of arrival can be measured to an accuracy of 72 ms in contrast to 
a value of 20 ms expected from the covariance matrix. As is well known, with the inclusion 
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of the post-Newtonian terms error in the estimation of the time of arrival and Newtonian 
chirp time increases by about a factor of 2 and 3, respectively |^ , ^ . 

As in the Newtonian case here again we see that the Monte Carlo curves approach the 
corresponding covariance matrix curves at a high SNR the only difference being that the 
agreement is reached at a higher SNR ~ 25. For SNRs larger than this the two curves 
are in perfect agreement with each other. As mentioned earlier, at^ shows an anomalous 
behaviour possibly arising out of insufficient resolution in the time of arrival and the chirp 
times. 

In Table ^ we have listed the true parameters A, the estimated values eA and the Monte 
Carlo errors ax for different SNRs. As in the Newtonian case here too the estimated values 
show a larger departure, than expected, from the true values. Histograms of the various 
measured parameters including tc are shown in Fig^ @ for a signal strength of 10. The 
skewness is below its standard deviation of y/Tb/Nt ]33[ | indicating the Gaussian nature of 
the various distributions. Even in the case of a post-Newtonian signal at^ is so small that 
we only have three non-zero bins in n(tc)- 

We now turn attention to other, more general, issues arising out of the simulations. 

In Sec. Ill C we have argued, on the basis of the behaviour of the noise-free correlation 
function, that the effective dimensionality of the parameter space for the purpose of de- 
tection, even in the case of a post-Newtonian signal, is only one-dimensional. The results 
of our Monte Carlo simulation unambiguously show that this is indeed true even in the 
presence of noise. We investigated the two-dimensional histogram, that gives the number of 
occurrances of different templates in the lattice, in a particular simulation. The templates 
that 'click' are all aligned along the line tq -t- ti = const. In a total of 5000 realisations 
corresponding to this simulation there is only one instance when a filter outside this region 
clicks giving a probability of less than 10"'^ for a template outside this region to give a 
maximum. Consequently, it is only necessary to choose a single filter along this curve. 

The distribution of the maximum correlation Cmax(A,tA) obtained from different noise 
realisations needs a special mention since it has an inherent bias. In Fig. ^ we have shown the 
distribution of the maximum correlation taken from one of our simulations corresponding to 
an SNR of 10. Notice a slight shift of the distribution towards a higher value and this cannot 
be accommodated within the expected fluctuation in the mean. The measured value of the 
standard deviation aj\^ is 0.95. Since the number of simulations is 5000 we expect that the 
signal strength should differ from the true value of 10 by no more than CT^/V5000 = 0.014. 
However, the mean value is 10.26 giving a deviation of 0.26 which is about 20 times larger 
than that expected. This occurs at all SNRs and for both Newtonian and post-Newtonian 
signals. This of course does not mean there is bug in the way we are computing the maximum 
correlation. In the process of maximisation, values greater than the signal strength are 
favoured and consequently the mean of the maximum correlation shows a shift towards a 
higher value. This suggests that that the maximum of the correlation is a biased estimator 
of the signal strength. We find, consistent with the covariance matrix calculation, that 
the amplitude parameter is uncorrelated with the rest of the parameters; crosscorrealtion 
coefficients DoM; M 0; (cf- 4-2) inferred from our Monte Carlo simulations are less 
than ~ 10"^. 

Finally, it is of interest to note how the phase parameter $ is correlated with the time 
of arrival. A plot of versus is shown in Fig. |l^. We find that the measured values 
of the time of arrival and the phase are such that 27r/o m^a = m*I^, where /o has a value of 
approximately 51 Hz. When the time of arrival shifts by more than a cycle of the signal 
the phase jumps by a factor 27r leading to the points seen in the top-left and bottom-right 
corner of the figure. This makes the estimation of the phase and the error in its estimation 
pretty involved. 



3. Incorporating the effects of upper cutoff 

As mentioned before, incorporating an upper cutoff at the onset of the plunge has a 
drastic effect on the estimation of parameters. The incorporation of the upper cutoff is 
implemented by stopping the waveform when the instantaneous frequency reaches the fre- 
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quency associated with the onset of the plunge or 750 Hz, whichever is lower. However, 
due to computational constraints we have carried out the simulations only for high mass 
binaries and hence the upper cutoff plays an important role in all our simulations. It is to 
be noted that the discussion of the ambiguity function in section HI C is not valid when the 
effect of the upper cutoff is taken into account though for low mass binaries, such as ns-ns 
binaries, the results there will still be valid. The further dependence of the signal waveform 
on the total mass of the system through the upper cutoff means that we can estimate the 
individual masses more exactly though the computational power is bound to increase. 

In order to carry out the simulations for the present case we selected a lOM© — 1.2Mq 
binary system as this enables us to choose the filter grid well within valid limits of tq and ri . 
The simulations were carried out for various values of SNR starting from 10. The histograms 
of the estimated parameters at an SNR of 10 is shown in fig. |ri|. At this SNR the errors 
obtained are arc — 39.3 ms, — 22.4 ms, at^ — 23.1 ms and at^ = 0.6 ms. These can be 
compared with the values in table ^ and we can see that except for the parameter tc the 
errors are substantially lesser when the upper cutoff is incorporated into the waveform. It is 
necessary to recompute the covariance matrix, as emphasized before, including the effect of 
the upper cutoff in order to compare these numerically obtained values with the covariance 
matri x. In order to do this it is not enough to replace the upper limit in the integral in eq. 
( |3.19| ) with the upper cutoff. The waveform now depends on the total mass of the system 
through the upper cutoff and this information has to be incorporated into the waveform. 

We carried out the simulations for various SNRs for the same value of masses quoted 
above. In the absence of the estimates of the covariance matrix when the upper cutoff is 
incorporated we assume that at an SNR of 40 the Monte Carlo estimates are consistent with 
those of the covariance matrix. In Figure we illustrate the results of our simulations. The 
points on the dotted line are the values obtained through Monte Carlo estimates and the 
continuous line is obtained by fitting a 1/p dependence of the errors on the SNR assuming 
consistency at an SNR of 40. It is seen as in the previous simulations that except for the 
parameter tc the other parameters are fairly consistent with a 1/p dependence of the errors 
on the SNR at an SNR greater than 25. 



V. CONCLUSIONS 



In this paper we have introduced the use of differential geometry in studying signal analysis 
and have addressed issues pertaining to optimal detection strategies of the chirp waveform. 
We have also carried out Monte Carlo simulations to check how well the covariance matrix 
estimates the errors in the parameters of the chirp waveform. We summarize below our 
main results. 

1. We have developed the concept of a signal manifold as a subset of a finite dimensional 
vector space of detector outputs. Using the correlation between two signal vectors 
as a scalar product we have induced a metric on the signal manifold. With this 
geometric picture it is possible to pose the question of optimal detection in a more 
general setting. We suggest that the set of template waveforms for the detection of the 
chirp signal need not correspond to any point on the chirp manifold. We propose an 
algorithm to choose templates off the signal manifold and show that the drop in the 
correlation due to the discreteness of the set of templates is reduced. This algorithm, 
though certainly not the best, motivates the search of more efficient templates. In 
addition, the chirp manifold corresponding to the second post-Newtonian waveform 
is shown to be effectively one-dimensional. This has important implications for the 
computational requirement for the on-line detection of the chirp signal. The use 
of a convenient set of parameters of the chirp waveform for carrying out numerical 
and analytical simulations is stressed. These parameters are such that the metric 
components are independent of the parameters which implies that the manifold is 
flat and the corresponding 'coordinate system' is Cartesian. As the metric defined is 
nothing but the Fisher information matrix, the covariance matrix, being the inverse 
of the Fisher information matrix, is also independent of the parameters. 
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2. Monte Carlo simulations have been carried out for the case of the initial LIGO to 
find out whether the actual errors in the estimation of parameters is consistent with 
the values predicted by the covariance matrix. Simulations have been carried out for 
both the Newtonian as well as the post-Newtonian waveforms. We have restricted 
ourselves to the case of high mass binary systems, such as bh-ns binaries, where the 
computational requirement is not very heavy since the length of the data train, in 
such cases, works out to be less than 8 s. Nevertheless, as has been shown in this 
paper, the covariance matrix is independent of the parameters identified by us when 
waveforms are terminated at a constant upper cutoff irrespective of their masses. 
Consequently, our results will hold good for binary systems of arbitrary masses. We 
point out the major problems that arise while performing a numerical simulation and, 
where appropriate, we suggest how they may be taken care of. In particular, the 
effect of incorporating the upper cutoff in the frequency of the gravitational wave at 
the onset of the plunge, which essentially depends on the total mass of the binary, 
is extremely important for high mass binaries. Since the covariance matrix with the 
inclusion of such a mass-dependent upper cutoff is not available we have carried out 
most of our simulations using a constant upper cutoff. This enables us to directly 
compare the results of our Monte Carlo simulations with those of the analytically 
computed covariance matrix. Since for binaries with total mass less than 5Mq the 
plunge induced upper cutoff is larger than that induced by the detector noise these 
effects can be ignored for such binaries. 

The numerical experiments indicate that the covariance matrix underestimates the 
errors in the determination of the parameters even at SNRs as high as 20. In the 
Newtonian case the correlation coefficient of the time of arrival ta and the Newtonian 
chirp time tq is found to be very close to —I, so much so that even at a SNR of 7.5, 
the instant of coalescence tc — ta + tq remains practically a constant. The error in 
the estimation of tq for the post-Newtonian waveform is about four times the error 
obtained in the case of the Newtonian waveform at the same SNR. This is expected 
as the first post-Newtonian correction to the waveform introduces a new parameter ti 
(called the first post-Newtonian chirp time) which is highly (anti)correlated with tq. 

For the post-Newtonian waveform at an SNR of 10 the error in tq is about 3 times 
that predicted by the covariance matrix. This corresponds to a factor of 2 in the 
chirp mass M. The distributions for the parameters have been obtained and arc seen 
to be unimodal distributions and are slightly more sharper than a Gaussian. When 
the plunge induced upper cutoff is incorporated into the waveform the errors in the 
estimation of parameters decrease by a factor of about 2.5. The correlation coefficient 
between tq and ti is also found to decrease, which is consistent with our discussion in 
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The results obtained suggest that higher moments than that used in obtained in com- 
puting the covariance matrix may be important in the determination of the errors in 
parameter estimation. In the geometric picture this amounts to taking into account 
curvature effects, either intrinsic or extrinsic. 

3. We suggest that tc is a more suitable parameter to estimate the direction to the source 
than the time of arrival. The latter is a kinematical parameter that fixes the time at 
which the gravitational wave frequency reaches the lower cutoff of the detector while 
the parameter tc has the physical significance of being the instant of coalescence. At 
an SNR of 10 the error in ta is too large (20 ms) to deduce the direction to the source 
accurately, whereas the error in the parameter tc is less than 0.5 ms. This will further 
go down substantially for the advanced LIGO. A detailed analysis of the determination 
of the direction to the binary using delays in tc between different detectors is carried 
out in Bhawal and Dhurandhar (also see [p9[). 

We now suggest further work which needs to be done along the lines of this paper. A 
full understanding of the chirp signal manifold when higher post-Newtonian corrections are 
incorporated into the waveform is in order. This will help in the development of more efficient 
algorithms for the choice of templates in the detection problem and facilitate reduction in 



25 



computational time. The Monte Carlo simulations which we have carried our are for the case 
of a binary waveform correct up to first post-Newtonian order. Moreover, only circular orbits 
are considered. The effect of eccentricity is currently being investigated Performing 
simulations when higher post-Newtonian corrections are taken into account calls for an 
immense amount of computational time. Fortunately, matched filtering algorithm being 
amenable to parallelization js^, one could aim at using the massively parallel computers, 
that are now becoming available the world over, in performing such simulations. 
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TABLE I. The estimated value of the parameters and their errors eA {a\) for the Newtonian 
waveform. The actual values of the parameters taken were fo = 5558.0 ms and ta = 200.0 ms. 
Except for the parameter A all values are quoted in ms. 

7.5 p = 10.0 p ^ 12.5 p ^ 15.0 p = 20.0 

eTo (aro) 5557 A (29.7) 5557.3 (15.6) 5557.6 (9.0) 5557.7 (6.2) 5558.5 (3.9) 

eta {at J 200.6 (29.4) 200.7 (15.5) 200.4 (8.9) 200.3 (6.2) 199.5 (3.9) 

etc (ate) 5758.0 (0.30) 5758.0 (0.20) 5758.0 (0.10) 5758.0 (0.08) 5758 (0.02) 

eA (ffA) 7.696 (0.96) 10.12 (0.98) 12.582 (0.99) 15.067 (0.99) 20.05 (0.99) 



TABLE II. The estimated value of the parameters and the errors in their estimation cA {a\) for 
the post-Newtonian waveform. The actual values of the parameters taken were to = 5558.0 ms, 
n = 684.0 ms and ta = 300.0 ms Except for the parameter A all values are quoted in ms. 

p = 10.0 p = 15.0 p = 20.0 p = 25.0 p = 30.0 



eTo (fro) 

cTl (cTti) 

eta (fta) 

etc (o-tc) 

eA (aA) 



5554.9 (136.1) 

685.55 (65) 
301.49 (73.26) 
6.542 (0.54) 
10.25 (0.97) 



5555.9 (64.8) 

685.2 (32.6) 
300.88 (33.4) 
6.542 (0.30) 
15.15 (0.98) 



5557.2 (30.1) 

684.6 (16.4) 
300.28 (14.5) 
6.542 (0.17) 
20.11 (.99) 



5556.7 (19.3) 

684.8 (10.5) 
300.5 (9.4) 

6.542 (0.10) 
20.1 (0.99) 



5558.7 (13.6) 

683.6 (7.2) 

299.7 (6.8) 
6.542 (0.06) 
30.1 (0.99) 
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FIG. 1. This figure illustrates the correlation function along the line of curvature as a function 
of To, the Newtonian chirp time, for the following three cases: (i) when the filter is placed on the 

manifold (dotted line), (ii) when the filter is the 'average' signal vector over the region (dashed 
line) and (iii) when the filter is chosen to be an appropriate linear combination of the previous two 
filters (solid line). 



FIG. 2. Contour diagram of the ambiguity function for the second post-Newtonian case. 

FIG. 3. Dependence of the errors in the estimation of the parameters of the post-Newtonian 
waveform on the upper cutoff frequency. The SNR is kept fixed at a value of 10. 

FIG. 4. Dependence of the errors in the estimation of the parameters of the post-Newtonian 
waveform on the total mass of the binary keeping the distance to the binary fixed. The waveforms 
are cutoff at frequencies corresponding to the onset of the plunge orbit and the SNR is normalized 
at a value of 10 for a 1.4 Mq-10 Mq binary. 

FIG. 5. Dependence of the errors in the estimation of parameters of the Newtonian waveform 
i.e. {o"ro , crt„ , crt(-, } as a function of SNR. The solid line represents the analytically computed errors 
whereas the dotted line represents the errors obtained through Monte Carlo simulations. 

FIG. 6. Distributions of the measured values of the parameters for the case of Newtonian signal. 
The total number of noise realizations is 5000. 

FIG. 7. Dependence of the errors in the estimation of parameters of the post-Newtonian waveform 

i.e. {cttq , (Tti , (Ttj, , CTtf, } as a function of SNR. The solid line represents the analytically computed 
errors whereas the dotted line represents the errors obtained through Monte Carlo simulations. 

FIG. 8. Distributions of the measured parameters of the post-Newtonian waveform at a SNR of 
25. The number of noise realizations is 5000. 

FIG. 9. The histogram of random variable mA at an SNR of 20. The variance in this parameter 
is independent of the SNR and is approximately equal to unity. 

FIG. 10. The correlation between ta and 3> is illustrated. The phase parameter simply follows 
the time of arrival parameter. 

FIG. 11. Distributions for the mcEisured parameters of the post- Newtonian waveform incorpo- 
rating the effect of the upper cutoff. The errors in the determination of the parameters is much 
smaller in this case. The total number of noise realizations is 5000. 

FIG. 12. The dependence of the errors on the SNR (p) when the upper cutoff is incorporated 
into the waveform 
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